R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[R.app GUI 1.52 (6188) i386-apple-darwin9.8.0]

[History restored from /Users/sr320/.Rapp.history]

> library( "DESeq" )
Loading required package: Biobase
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following object(s) are masked from ‘package:stats’:

    xtabs

The following object(s) are masked from ‘package:base’:

    anyDuplicated, cbind, colnames, duplicated, eval, Filter,
    Find, get, intersect, lapply, Map, mapply, mget, order,
    paste, pmax, pmax.int, pmin, pmin.int, Position, rbind,
    Reduce, rep.int, rownames, sapply, setdiff, table, tapply,
    union, unique

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages
    'citation("pkgname")'.

Loading required package: locfit
locfit 1.5-8 	 2012-04-25
> 
> 
> CountData<- read.table("Gill_GonadM_v9.txt",row.names=1, header=TRUE,)
> 
> head (CountData)
             Gill Gonad_M
CGI_10000001    0       0
CGI_10000002    7       6
CGI_10000003    0       0
CGI_10000004    0       0
CGI_10000005    0       0
CGI_10000009   65      50
> 
> Treatment <- factor( c("Gill","Gonad_M") )
> 
> cds <- newCountDataSet( CountData, Treatment )
> 
> LibrarySize <- estimateSizeFactors( cds )
> sizeFactors( LibrarySize )
     Gill   Gonad_M 
0.8363982 1.1956027 
> Dispersons <- estimateDispersions( LibrarySize, method="blind", sharingMode="fit-only" )
Warning message:
In parametricDispersionFit(means, disps) : Dispersion fit did not converge.
> DE <- nbinomTest( Dispersons, "Lean", "Siscowet" )
Error in if (dispTable(cds)[condA] == "blind" || dispTable(cds)[condB] ==  : 
  missing value where TRUE/FALSE needed
> 
> DE <- nbinomTest( Dispersons, "Gill", "Gonad_M" )
> 
> head (DE)
            id  baseMean baseMeanA baseMeanB foldChange log2FoldChange
1 CGI_10000001  0.000000  0.000000  0.000000        NaN            NaN
2 CGI_10000002  6.693804  8.369219  5.018389  0.5996246     -0.7378686
3 CGI_10000003  0.000000  0.000000  0.000000        NaN            NaN
4 CGI_10000004  0.000000  0.000000  0.000000        NaN            NaN
5 CGI_10000005  0.000000  0.000000  0.000000        NaN            NaN
6 CGI_10000009 59.767044 77.714178 41.819911  0.5381246     -0.8939878
       pval padj
1        NA   NA
2 0.9616185    1
3        NA   NA
4        NA   NA
5        NA   NA
6 0.7714456    1
> 
> plotDE <- function (DE)
+ 	plot(
+ 		DE$baseMean, 
+ 		DE$log2FoldChange,
+ 		log="x", pch=20, cex=.3,
+ 		col = ifelse (DE$padj < .05, "red", "black"))
> 		
> plotDE(DE)	
Warning message:
In xy.coords(x, y, xlabel, ylabel, log) :
  4634 x values <= 0 omitted from logarithmic plot
> hist(DE$pval, breaks=100, col="skyblue", border="slateblue", main="")
> 
> write.table(DE, "Oyster_GillGondad_NATURE_DESeq.txt", row.names = FALSE, sep="\t")
> 
>